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ABSTRACT 


We  have  constructed  a  global-scale  model  of  P-wave  velocity  with  an  emphasis  on  improving  travel  time 
prediction  at  both  regional  and  teleseismic  distances  simultaneously.  The  LLNL-G3Dv2  tomographic  model 
is  built  within  a  spherical  tessellation  framework  whereby  irregular  and  discontinuous  surfaces  are 
explicitly  represented.  Fully  3-D  ray  tracing  is  employed  for  travel  time  prediction.  The  data  consist  of 
-2.7  million  P  and  Pn  arrivals  that  are  re-processed  using  our  global  multi-event  locator  known  as 
Bayesloc.  Bayesloc  is  a  formulation  of  the  joint  probability  distribution  across  multiple-event  location 
parameters,  including  hypocenters,  travel  time  corrections,  pick  precision,  and  phase  labels.  Modeling  the 
whole  multiple-event  system  results  in  accurate  locations  and  an  internally  consistent  data  set  that  is  ideal 
for  tomography.  Our  recently  developed  inversion  approach  called  Progressive  Multi-level  Tessellation 
Inversion  (PMTI)  captures  regional  trends  and  fine  details  where  data  warrant.  Using  PMTI,  we  model 
multiple  heterogeneity  scale  lengths  without  resorting  to  user-defined  parameter  grids  that  can  affect  the 
resulting  model.  A  number  of  features  have  emerged  from  the  images  including  details  of  the  underthrusted 
Arabian  lithosphere  beneath  Eurasia  (shown  in  our  previous  reports)  and  sharp  details  of  subducted  slabs 
around  the  world.  Based  on  preliminary  relocation  tests,  using  LLNL-G3Dv2  travel-time  predictions 
typically  reduces  mislocation  error  for  a  set  of  globally  distributed  GT0-5  events  by  30-45%  (from  -10- 
1 1km  using  akl  3 5  to  -6-7  km)  when  a  sufficient  number  of  data  are  used. 
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OBJECTIVES 


The  objective  of  this  project  is  to  generate  a  seamless  model  of  the  Earth’s  crust  and  mantle  that  is  capable 
of  accurately  predicting  regional  and  teleseismic  travel  times.  The  objective  will  further  advance  event 
monitoring  capabilities,  through  self-consistent  determination  of  seismic  travel  times  that  are  necessary  for 
accurate  event  location.  The  objective  requires  simultaneous  modeling  of  detailed  upper  mantle 
heterogeneities  to  adequately  predict  regional  phases  such  as  Pn  and  lower  mantle  structures  to  predict 
teleseismic  travel  times.  In  addition  to  determining  velocity  structure,  an  efficient  model  design  and 
computational  framework  must  be  established  in  order  to  routinely  utilize  the  outcome  model.  In  this 
report,  we  present  our  latest  global  3-D  P-wave  velocity  model  and  event  location  validation  results. 

RESEARCH  ACCOMPLISHED 


Employing  a  number  of  custom  procedures  described  in  previous  Monitoring  Research  Review 
proceedings  (e.g.,  Simmons  et  al.  2009a,  2010a)  and  the  most  recent  peer-reviewed  publications  (e.g., 
Myers  et  al.  2011;  Simmons  et  al.  2011),  we  have  developed  a  revised  global-scale  P-wave  velocity  model 
that  incorporates  -2.7  million  re-processed  regional  and  teleseismic  travel  times.  We  also  performed 
exhaustive  event  location  validation  tests  for  a  globally  distributed  set  of  events  with  Ground  Truth  (GT) 
levels  of  5  or  less.  With  this  validation  set  of  events,  we  typically  find  30-45%  improvement  in  median 
epicenter  error  when  determining  location  on  the  basis  of  the  most  recent  3-D  model  ( LLNL-G3Dv2 )  when 
more  than  -20  arrivals  are  available.  Although  typical,  this  level  of  improvement  depends  on  the  number  of 
supporting  events  and  stations  in  a  particular  region  (i.e.  tomographic  resolution). 


Figure  1.  The  dataset  consists  of  -13,400  seismic  events  (in  blue)  and  -8,000  stations  (in  green) 

providing  a  total  of  3.4  million  arrivals  for  phases:  P,  Pn,  PcP,  pP,  sP,  Sn,  Pg,  and  Lg.  The 
data  are  culled  by  the  Bayesloc  process  and  -2.7  million  P  and  Pn  phases  are  used  to 
develop  the  tomography  model  ( LLNL-G3Dv2 ). 

A  New  Global  Dataset 

Our  raw  (initial)  compilation  of  travel  time  measurements  come  from  a  subset  of  our  local  database  at 
Lawrence  Livermore  National  Laboratory  (Ruppert  et  al.,  2005)  and  the  publicly  available  Engdahl-van  der 
Hilst-Buland  (EHB)  catalog  (Engdahl  et  al.,  1998).  The  most  well-recorded  global  seismic  events  were 
selected  nominally  within  1°  lateral  bins  and  within  6  depth  bins  at  35,  75,  150,  300,  450  and  700  km 
depths.  The  dataset  includes  all  GT5  and  better  events  with  >10  teleseismic  or  regional  arrival  picks.  In 
total,  there  are  approximately  3.4  million  arrivals  from  -13,400  seismic  events  recorded  at  -8,000  stations 
around  the  globe  (figure  1).  The  vast  majority  of  the  data  are  teleseismic  P-wave  arrivals,  but  we  also 
include  the  phases:  Pn,  PcP,  pP,  sP,  Sn,  Pg  and  Lg. 
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We  have  extended  the  Bayesloc  method  (originally  designed  to  simultaneously  locate  a  cluster  of  events)  to 
global-scale  multiple-event  relocation  (see  Myers  et  al.  2007,  2009,  2011).  The  Bayesloc  method 
simultaneously  models  all  event  hypocenters,  regional  travel  time  biases,  path  specific  travel  time  errors, 
arrival  time  measurements  and  arrival  time  phase  assignments  within  a  Bayesian  statistical  framework.  The 
process  results  in  accurate  event  locations  and  travel  times  which  are  essential  when  generating  a  travel 
time  tomography  model  (see  Simmons  et  al.,  2011).  The  updated  Bayesloc  procedure  was  performed  on  the 
aforementioned  set  of  13,400  events  and  3.4  million  arrivals  to  produce  a  new  set  of  data  for  tomography 
with  revised  event  locations  (Figure  2).  In  summary,  we  find  a  median  epicenter  shift  of  25.9  km  from  the 
bulletin  locations.  After  Bayesloc  processing,  the  travel  time  residual  standard  deviation  drops  from  1.87  to 
1.35  seconds  with  respect  to  the  AK135  model  (Figure  3).  Perhaps  most  importantly,  we  find  clear  regional 
trends  in  epicenter  shift  which  has  substantial  impact  on  the  resulting  tomographic  model. 


Figure  2.  Epicenter  shift  vectors  demonstrating  the  lateral  re-location  differences  from  the  raw  data 
and  the  Bayesloc  processed  data  set.  The  Bayesloc  process  results  in  a  median  epicenter 
shift  of  25.9  km  with  clear  regional  trends. 
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Figure  3.  Travel  time  residual  distribution  for  the  ~2.7  million  P  and  Pn  arrivals.  (Left)  The  initial 
bulletin  data  with  respect  to  the  AK135  model.  (Center)  The  Bayesloc  processed  travel 
times  with  respect  to  the  akl35  model.  (Right)  The  Bayesloc  processed  travel  times  with 
respect  to  the  LLNL-G3Dv2  model  presented  in  this  report. 


The  LLNL-G3Dv2  Model 

The  LLNL-G3D  series  of  models  are  designed  with  the  use  of  spherical  tessellations  whereby  irregular  and 
discontinuous  surfaces  can  be  explicitly  represented  (Ballard  et  al.,  2009;  Simmons  et  al.,  2011).  This  style 
of  explicit  Earth  representation  is  desirable  given  that  we  are  simultaneously  modeling  regional-  and 
global-scale  velocity  structures  within  a  single,  self-consistent  model.  Specifically,  3-D  ray  tracing  is 
performed  due  to  the  extreme  path- velocity  dependencies  arising  at  regional  distances.  Thus,  we  determine 
sensitivity  via  3-D  ray  tracing  (based  on  Zhao  et  al.,  1992)  while  considering  irregular  discontinuity 
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structure  and  continuous  media  simultaneously.  In  addition,  sensitivity  is  spread  across  broad  depth  zones 
and/or  multiple  model  units  to  mitigate  the  issue  of  path-velocity  interdependence  as  well  as  severe  multi- 
pathing.  See  Simmons  et  al.  (2011)  for  a  complete  description  of  the  3-D  ray  tracing  procedures  employed 
within  the  current  study. 

Rather  than  using  a  simple  1  -D  velocity  model,  we  chose  to  create  a  starting  P  wave  velocity  model  by 
combining  components  of  existing  models  of  the  crust  and  mantle.  There  are  many  benefits  of  leveraging 
past  work  to  design  a  first  order  starting  model  including:  1)  providing  constraints  on  regions  with  poor  P- 
wave  coverage,  2)  mitigation  of  the  non-linear  relationship  between  ray  paths  and  velocity  structure,  and 
3)  allowing  for  the  determination  of  a  new  model  that  is  most  consistent  with  what  is  generally  understood 
about  the  3-D  structure  of  the  crust  and  mantle  from  decades  of  research. 

For  the  crust,  we  use  the  ‘Unified’  crust  model  from  the  work  of  Pasyanos  et  al.  (2004),  Steck  et  al.  (2004), 
and  Bassin  et  al.  (2000).  The  crust  model  consists  of  7  units  including  water,  3  sediment  layers,  and  3 
crystalline  crust  layers.  For  the  shallowest  upper  mantle,  we  use  the  results  of  the  Regional  Seismic  Travel 
Time  (RSTT)  model  described  in  Myers  et  al.  (2010)  including  the  sub-crustal  velocity  and  mantle  gradient 
terms.  For  the  remainder  of  the  mantle,  we  incorporated  the  results  of  the  GyPSuM  model  (Simmons  et  al. 
2010b).  GyPSuM  is  a  3-D  model  of  mantle  S  wave  speed,  P  wave  speed  and  density  developed  through 
simultaneous  inversion  of  seismic  and  a  suite  of  geodynamic  constraints.  The  GyPSuM  model  represents 
the  latest  of  a  sequence  of  joint  global  modeling  efforts  described  in  Simmons  et  al.  (2006,  2007,  2009b). 

Numerous  studies  have  independently  concluded  that  significant  topographic  variations  of  the  transition 
zone  discontinuities  exist.  More  specifically,  it  has  been  shown  that  the  410  km  and  660  km  discontinuities 
vary  in  depth  by  ±30  km  or  more  over  relatively  short  length  scales  (Lawrence  and  Shearer  2008).  These 
variations  may  lead  to  incorrect  3-D  ray  path  geometries,  incorrect  model  sensitivity  estimates,  and 
therefore  inaccurate  velocity  structure  after  tomographic  inversion  is  performed.  For  this  reason,  we 
adjusted  the  depths  of  the  410  and  660  according  to  the  transition  zone  discontinuity  topographies 
determined  in  the  global  high-resolution  SS  precursor  study  of  Lawrence  and  Shearer  (2008).  In  addition, 
the  model  explicitly  characterizes  Earth’s  asphericity  in  accord  with  the  WGS84  reference  ellipsoid  and  the 
expected  hydrostatic  shape  of  mantle  layers  and  the  core-mantle  boundary. 


Juan  de  Fuca  /  Gorda  Cocos 

.  •••«  .••• 


Figure  4.  Cross-sections  through  the  LLNL-G3Dv2  model  showing  subducting  slabs  beneath  the 
Americas. 
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As  noted  in  previous  reports,  the  tessellation-based  model  design  is  hierarchical  since  the  mesh  is  a  product 
of  subdividing  a  base  level  tessellation  a  number  of  times.  One  of  the  major  outstanding  issues  in  seismic 
tomography  is  the  uneven  sampling  of  seismic  data  and  the  differing  wavelengths  of  actual  seismic 
heterogeneity.  These  issues  make  it  difficult  to  design  an  appropriate  mesh.  Thus,  we  have  developed  an 
inversion  process  (PMTI)  that  exploits  the  hierarchical  nature  of  the  tessellation-based  design  and  allows 
the  data  to  determine  the  level  of  model  complexity,  without  user  bias.  PMTI  serves  as  an  alternative  to 
existing  multiresolution  approaches  and  we  have  demonstrated  that  the  process  robustly  images  regional 
trends  while  allowing  localized  details  to  emerge  where  resolution  is  sufficient  (see  Simmons  et  al.  2011 
for  more  details).  Using  the  PMTI  process,  we  update  the  aforementioned  starting  model  to  produce  the 
P-wave  velocity  model  shown  in  Figures  4-5. 


Philippines 


Marianas 


Figure  5.  Cross-sections  through  the  LLNL-G3Dv2  model  showing  subducting  slabs  beneath  Eurasia 
and  Indonesia,  as  well  as  the  superplume  beneath  Africa.  Color  scales  change  from  one 
frame  to  the  next  (denoted  by  the  colorbars)  and  all  values  are  in  terms  of  percent  of 
compressional  velocities  (Vp)  variation  relative  to  the  mean  model  with  depth. 
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The  resulting  image  (Figure  4-5)  reveals  a  number  of  subducted  slabs  around  the  world,  including:  i)  details 
of  the  Juan  de  Fuca,  Cocos,  and  Nazca  plates  beneath  the  Americas,  ii)  the  Pacific  and  Phillipine  plates 
along  the  western  margins  of  the  Pacific,  iii)  the  African  plate  subducting  beneath  Europe  along  the 
Hellenic  Arc,  iv)  interaction  of  the  Australian,  Phillipine  and  Pacific  plates  beneath  Indonesia,  and 
v)  superplumes  beneath  the  Pacific  and  Africa.  Countless  other  features  emerge  from  the  images  but  further 
discussion  is  beyond  the  scope  of  this  report. 

Resolution  tests  were  carried  out  using  the  inversion  procedures  discussed  above.  Our  synthetic  model  is  a 
combination  of  two  checkerboards  summed  together  to  create  large-scale  regional  trends  combined  with 
more  local  detail  (Figure  6).  The  smallest  checkerboard  squares  are  5°  x  5°  in  dimension  and  the  pattern 
was  repeated  with  opposite  signs  in  each  of  the  3 1  inversion  layers  from  the  Moho  to  the  core-mantle 
boundary.  Therefore,  we  developed  a  single  model  for  the  entire  depth  range  of  the  mantle  rather  than 
testing  one  layer  at  a  time  (Figure  6  and  the  left  column  of  Figure  7).  In  addition,  we  carried  out  resolution 
tests  with  a  single  layer  at  a  time  approach  which  is  the  standard  approach  to  compare  the  two  styles 
(Figure  7).  Not  surprisingly,  the  “one  layer  at  a  time”  approach  shows  better  input  model  recovery 
especially  where  data  densities  are  low  (e.g  beneath  the  Pacific  Ocean  and  near  the  South  Pole). 


Figure  6.  Checkerboard  test  input  pattern,  recovery  at  the  Moho,  and  recovery  along  the  green,  360- 
degree  cross-section. 
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Multiple  Layers  Simultaneously 
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Figure  7.  Resolution  test  recovery  at  3  specific  depths  for  the  case  with  multiple  layers  (left  column) 
and  the  case  with  a  single  layer  at  a  time  (right  column)  for  comparison. 


Global  Event  Location 

For  location  validation,  we  selected  1 16  of  the  most  well-recorded  GT5  (or  better)  events  within  5-degree 
bins.  Each  event  that  maximizes  the  number  of  arrivals  or  network  coverage  (minimizes  station  azimuthal 
gap)  at  regional  or  teleseismic  distances  was  selected.  As  a  result,  up  to  4  distinct  events  may  be  selected 
within  each  5 -degree  bin.  The  events  are  globally  distributed,  but  the  majority  of  events  are  in  Eurasia 
(Figure  8).  Data  from  the  116  events  were  excluded  from  the  tomographic  inversion  to  produce  a  validation 
version  of  LLNL-G3Dv2  to  prevent  circularity.  For  each  event,  we  randomly  selected  specific  numbers  of  P 
and  Pn  arrival  times  to  be  used  for  event  location  (ranging  from  0  to  100  arrivals).  For  each  event  data 
sampling,  we  required  that  there  be  at  least  5  arrivals  and  that  enough  data  were  available  to  produce  at 
least  10  distinct  realizations  of  arrival  times. 


All  relocations  were  determined  using  the  Bayesloc  algorithm  in  “single-event”  mode  -  no  travel  time 
corrections,  assessment  of  travel  time  error,  or  phase  labeling  error  -  and  without  significant  assumed  prior 
information  about  the  event.  For  the  1-D  model  case,  locations  were  determined  based  on  akl  3 5  travel 
times.  For  the  3-D  case,  we  used  an  iterative  approach  with  the  steps:  1)  find  the  1-D  location,  2)  adjust  the 
arrival  times  based  on  the  difference  between  akl 3 5  and  the  computed  LLNL-G3Dv2  travel  times, 
and  3)  relocate  using  Bayesloc.  Effectively,  this  process  removes  the  3-D  travel  time  signal  from  the  data 
given  a  location  that  is  “in  the  ballpark”  (i.e.,  1-D  location)  and  then  relocate  once  again  assuming  a 
1-D  background  velocity  model  (akl  3 5). 
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Figure  9  shows  a  summary  of  the  event  location  validation  results.  The  results  are  demonstrated  by  binning 
the  event  realizations  according  to  the  number  of  P  and  Pn  used,  and  the  median  epicenter  mislocation  error 
(in  km)  was  determined  for  all  events  located  with  the  specified  number  of  data.  We  find  that,  in  most 
cases,  using  the  3-D  model  produces  more  accurate  locations  than  the  akl 3 5  model.  More  specifically,  in 
98%  of  all  the  cases,  the  3-D  model  improved  location  accuracy.  The  exceptions  occur  for  relocations  using 
few  data,  where  both  1-D  and  3D  models  result  in  mislocations  of  many  tens  of  km. 

Aside  from  cases  with  very  few  data,  we  find  the  least  improvement  when  only  P  arrivals  are  used 
(Figure  10).  With  sufficient  data,  the  median  mislocation  using  akl 3 5  is  -9-10  km  whereas  LLNL-G3D 
mislocates  by  -7  km.  The  -20-30%  improvement  is  not  surprising  since  akl  3 5  is  a  good  global  average 
model  for  teleseismic  arrivals  by  design.  In  cases  where  only  Pn  arrivals  are  used,  the  3-D  model  improves 
location  more  considerably  more  (-11  km  to  -6-7  km,  -35-45%)  owing  to  significant  3-D  heterogeneities 
on  the  regional  scale. 


Figure  8.  Events  used  for  location  validation.  See  text  for  details. 


Figure  9.  Median  epicenter  error  for  all  of  the  realized  validation  events  (see  text).  (Left)  Epicenter 
errors  as  a  function  of  the  number  of  P  and  Pn  picks  used  assuming  the  AK135  1-D  velocity 
model.  (Right)  Using  the  LLNL-G3Dv2  model. 
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Number  of  P  picks 


Number  of  Pn  picks 


Figure  10.  Median  epicenter  mislocation  for  the  cases  where  all  data  are  either  P  picks  (left)  or  Pn 
picks  (right)  where  the  other  phase  is  excluded. 


CONCLUSIONS  AND  RECOMMENDATIONS 


We  have  constructed  a  new  version  of  the  LLNL-G3D  (version  2)  global-scale  model  of  P-wave  velocity 
that  is  a  seamless  model  of  Earth’s  crust  and  upper  mantle.  The  LLNL-G3Dv2  tomographic  model  is 
constructed  with  -2.7  million  P  and  Pn  arrivals  that  are  re-processed  using  our  global  multi-event  locator 
known  as  Bayesloc.  With  LLNL-G3Dv2 ,  the  large  set  of  globally  distributed  P  wave  data  are  predicted 
within  1.08  seconds  standard  deviation  compared  to  1.87  seconds  for  akl35  (-67%  variance  reduction). 

The  image  provides  sharp  details  of  subducted  slabs  around  the  world,  including  those  most  often  only  seen 
with  detailed  regional  studies,  as  well  superplumes  and  many  other  geologically  significant  features. 

Based  on  these  preliminary  tests,  LLNL-G3Dv2  typically  reduces  mislocation  error  for  a  set  of  1 16  globally 
distributed  GT0-5  events  by  30-45%  (from  -10-1 1km  using  akl35  to  -6-7  km)  when  a  sufficient  number 
of  data  are  available.  The  3-D  model  predicts  locations  that  are  more  accurate  than  akl 3 5  in  98%  of  the 
cases  considered.  Most  notably,  the  3-D  model  more  substantially  improves  the  location  of  events  when  no 
P  wave  arrivals  are  used  in  the  location  (Pn  only).  This  result  is  expected  due  to  the  fact  that  akl 3 5  was 
designed  to  be  an  average  Earth  model  to  compute  average  P-wave  velocity  times  at  teleseismic  distances. 

We  find  that  regions  with  low  seismicity  levels  and  light  station  coverage  produce  the  worst  location 
predictions  with  the  3-D  model  owing  to  the  limited  understanding  of  the  velocity  structures.  The  extreme 
example  is  in  Northwest  Asia/Eastern  Europe  (typically  -18%  epicenter  location  improvement)  compared 
to  other  regions  such  as  western  North  America  (typically  >45%  location  improvement).  The  limitation  of  a 
specific  type  of  seismic  information  (e.g.,  P  and  Pn  arrivals)  in  a  particular  region  suggest  that  other  forms 
of  data  must  be  incorporated  when  generating  a  travel  time  model  suitable  for  accurate  event  monitoring. 
This  should  be  the  focus  of  future  work. 
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